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In modeling nonequilibrium systems one usually starts with a definition of the microscopic dy- 
namics, e.g., in terms of transition rates, and then derives the resulting macroscopic behavior. We 
address the inverse question for a class of steady state systems, namely complex fluids under con- 
tinuous shear flow: how does an externally imposed shear current affect the microscopic dynamics 
of the fluid? The answer can be formulated in the form of invariant quantities, exact relations for 
the transition rates in the nonequilibrium steady state, as discussed in a recent letter [A. Baule 
and R. M. L. Evans, Phys. Rev. Lett. 101, 240601 (2008)]. Here, we present a more pedagogical 
account of the invariant quantities and the theory underlying them, known as the nonequilibrium 
counterpart to detailed balance (NCDB). Furthermore, we investigate the relationship between the 
transition rates and the shear current in the steady state. We show that a fluctuation relation of 
the Gallavotti-Cohen type holds for systems satisfying NCDB. 

PACS numbers: 05.70.Ln, 05.20.Jj, 83.60.Rs 

I. INTRODUCTION 

One of the greatest triumphs of nineteenth century physics was the development of equilibrium statistical mechanics 
- a theory that exactly determines the probability of encountering any given arrangement of (the positions and 
momenta of) a large collection of interacting bodies (atoms, molecules, etc.), even though their equations of motion 
are intractable. From that probabilistic description, macroscopic properties en masse can be found, such as a fluid's 
pressure, density, anisotropy, etc. The theory assumes the existence of thermodynamic equilibrium, an asymptotic 
long-time limit solution to the entire dynamics in which the statistical properties are unchanging and independent 
of initial conditions. Conversely, since the birth of statistical mechanics an exact prediction of the statistics of 
the particles' motion has been unobtainable in the presence of flow, instead requiring various approximations or 
simplifications, which are the bread and butter of much present-day research into nonequilibrium phenomena. 

There is a certain class of highly nonequilibrium states of matter that possesses time-invariant statistics, and 
therefore might be approachable by exact statistical analysis without simplified modeling or near-equilibrium approx- 
imations. Such states arise when continuous shear flow is applied to complex fluids — materials composed of classical 
particles much larger than atoms. Under flow the structure of complex fluids can be radically re-ordered, even to the 
point of undergoing shear-induced phase transitions. 

For instance, a solution of amphiphiles in the worm-like micelle phase can exhibit a negative gradient in the 
characteristic flow curve if driven at high shear rates. This means that for a range of shear rates the viscosity decreases 
under shear. If the fluid is continuously sheared in this regime, it becomes mechanically unstable and separates into 
'bands' of high and low shear rate Also, at higher concentration the solution of amphiphiles can undergo a 
phase transition from a lamellar phase into complex densely packed onion-like structures [2J. Both phenomena 
are structural phase transitions that are controlled by shear rate in addition to the usual parameters (temperature, 
pressure, concentration) and are examples in which the fluid undergoes drastic microscopic re-configuration in response 
to the imposed (macroscopic) driving at the boundaries. 

A fluid in a sheared steady state is described by the same Hamiltonian as at equilibrium, since no external field 
is applied to drive the system; only the boundary conditions are different. Nevertheless, since they contain non-zero 
fluxes, such states fall outside the jurisdiction of equilibrium statistical mechanics, despite exhibiting all the vast 
variety of reproducible behaviors, structures, and transitions seen in equilibrium thermodynamic systems. 

A statistical description of complex fluids under continuous shear flow, starting from first principles, is provided by 
a theory of nonequilibrium transition rates, known as the nonequilibrium counterpart to detailed balance (NCDB) 
S EL 0- This theory addresses the following question: given an externally imposed macroscopic shear current at 
the boundaries of the fluid, what is the effect on the microscopic dynamics (the transition rates) in the bulk of the 
fluid? In the absence of any driving, that is, at thermal equilibrium, the microscopic transition rates are constrained 
to satisfy the principle of detailed balance, namely that the ratio of forward to reverse transition rates between any 
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pair of microstates must equal the Boltzmann factor of their energy difference: 

UJji 



(1) 



where w, 3 denotes the transition rate from microstate i to j and (3 the inverse temperature. Detailed balance is 
a consequence of the influence of the fluid's thermal surroundings: it is immersed in a larger volume of fluid at a 
certain temperature which acts as an equilibrium heat reservoir. The stochastic influence of this heat reservoir puts 
constraints on the allowed microscopic transitions in the form of Eq. (TT]) . 

If shear is imposed at the boundaries of the heat reservoir detailed balance no longer holds. However, a fluid region 
within the bulk continues to receive stochastic forces from the reservoir, which is now itself under flow. In the steady 
state this nonequilibrium heat reservoir imposes constraints on the transition rates according to NCDB, which can be 
expressed in the form of a one-to-one mapping between the transition rates at equilibrium and those in the sheared 



Here, f^j (y) denotes a transition rate in the sheared steady state, parametrized by v, a Lagrange multiplier character- 
izing the driving strength of the reservoir. Under the influence of shear forces applied at the boundaries of the whole 
ensemble, the rates are thus enhanced or attenuated with respect to equilibrium. The factor e" Axji exceeds unity if 
the transition i — > j involves a conformational change that increments the shear strain by a positive amount Axij. 
This factor simply boosts every transition in the forward direction irrespective of the state space structure. By itself, it 
would represent a simple mean-field expression for the driven transition rates. Important non-mean-field information 
about the global properties of the state space is contained in the quantities Aqji (y) , whose precise definition is given 
further below, in Sec. IIIIl 

Recently, it has been demonstrated that the constraints of NCDB can be cast into the form of invariant quantities 
that remain unchanged by the driving Q and apply to any pair of microstates in the following way: (i) The product of 
forward and reverse transition rates remains invariant under the driving, i.e., is the same in the equilibrium and in the 
sheared steady state: ujijUji = Qijttji (where fiy denotes the transition rate in the sheared state), (ii) The difference 
of total exit rates remains invariant: (Qjfc — Cljk) = X^fc {^ik ~~ u jk)~ These invariant quantities represent exact 
relations for the transition rates in the sheared steady state, arbitrarily far away from equilibrium and represent a 
prediction of NCDB that can be tested in an experiment. Furthermore, using a graph representation of NCDB we have 
devised a systematic method to calculate the driven transition rates from a set of algebraic equations for arbitrary 
configurations of the system state space, thus greatly enhancing the applicability of the theory. 

In this article we present a more pedagogical account of these recent results and the theory of NCDB underlying 
them. In particular, we provide a detailed derivation of the various representations of NCDB from a nonequilbrium 
sheared ensemble. Furthermore, we investigate the properties of the shear current in systems satisfying NCDB 
and show that the shear current exhibited by an individual system trajectory satisfies a fluctuation relation of the 
Gallavotti-Cohen type. 

The remainder of this article is organized as follows. The nonequilibrium ensemble on which the derivation of 
NCDB relies is introduced in Sec.[Hj In Sec. IIIIl we present a detailed review of the different representations of NCDB 
as previously discussed in Refs. [311]) III- A graph representation is presented in Sec. IIVI providing an intuitive way 
to discuss Master equation systems in discrete state spaces. The invariant quantities are derived in Sec. El where we 
also formulate the systematic calculation method for the transition rates in the sheared steady state. The relationship 
between the shear current and the transition rates is discussed in Sec. I VII Our rather formal results are elucidated 
in two simple hopping models which allow for an explicit calculation of the driven transition rates and other relevant 
quantities of the NCDB formalism (Sees. IVTlI and lyTIi]) . Even in these simple models it is evident that the predictive 
power of NCDB goes well beyond simple mean-field theories. 



In order to describe nonequilibrium states of complex systems subject to noise, one usually relies on a probabilistic 
description in terms of transition rates. In such an approach the main quantity of interest is the set of probability 
distributions {pi{t)} over states i = l,...,n, which express the probability to find the system in state i at time t. The 
nature of these states depends on the level of description; for reasons that will become clear below they are considered 
to be classical microstates in the following. The dynamical evolution of the system is then governed by the Master 
equation, a balance equation for the probability: 



steady state [1] 



(2) 



II. NONEQUILIBRIUM ENSEMBLE AND PATH ENTROPY 




(3) 
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The sum is here taken over the set of states {j} connected with i, where w, 3 denotes the rate of transition from state 
i to state j. The difference u>jiPj(t) — u>ijPi(t) is interpreted as the microscopic probability current between states i 
and j, and therefore Eq. |3]) states the conservation of probability. For a probabilistic interpretation of the Pi(t) we 
require the conditions < Pi(t) < 1 and J2i=iPi(t) = 1- 

A steady state is characterized by stationarity of the statistics, implying that all single-time distributions are time 
independent. If we set dpi(t)/dt — in Eq. ([3|) we find that the condition for a steady state is given by the balance 
of total in- and outgoing flow for every state i 



Nonequilibrium steady states are characterized by a non-zero net flow of particles, heat, etc., running through the 
system. By contrast, in equilibrium this flow vanishes, a fact which is manifest as the strong condition of detailed 
balance: the net probability current between any two configurations of the system is zero at equilibrium. This is 
expressed as 



for all i, j so that Eq. is trivially satisfied. 

The condition of detailed balance is a physical property of transition rates in equilibrium systems and can be derived 
from Liouville's equation for closed, isolated, physical systems [7]. 'Physical' refers to the existence of an underlying 
Hamiltonian such that the dynamics is deterministic and microscopically reversible. 'Isolation' implies that energy is 
a conserved quantity, i.e., a constant of motion. The trajectory in phase-space is then restricted to a shell of constant 
energy and no transitions between shells of different energy are possible. For a canonical system, these considerations 
apply in the same way to the combination system plus reservoir. The probability distribution of the system is then 
known to satisfy Boltzmann's law p^ q oc e~^ Ei . In this context detailed balance can be interpreted as a statement of 
four fundamental properties characterizing system and heat reservoir: (i) ergodicity, (ii) microscopic reversibility, (iii) 
time-translation invariance of statistical properties, and (iv) conservation of energy. As a result of these properties 
there are m/2 constraints acting on m transition rates of the system, expressed by Eq. (fTj). 

Having specified an equilibrium state in this way, we may consider a particular class of driven steady state systems, 
namely fluids under continuous shear. Let us consider a fluid region far from the boundaries as our 'system'. The heat 
reservoir consists of the fluid volume surrounding this region. If we assume that any correlation lengths are negligibly 
small compared with the fluid volume (a condition that may be unenforcable for turbulent flows), then this reservoir 
is only characterized by its macroscopic observables, which are in this case mean energy and mean shear rate. The 
stochastic influence on the system is here not that of an equilibrium reservoir, but of a reservoir which is itself in 
a nonequilibrium condition, under shear. In this case the properties (i) — (iii) remain valid. Ergodicity is generally 
difficult to prove rigorously even for the simplest systems at equilibrium, but it is assumed on empirical grounds, since 
experimental observations are repeatable, irrespective of precise initial conditions. Since the shear acts only at the 
boundaries of the reservoir, the dynamics of individual molecules is still governed by the same equations of motion as 
in equilibrium and thus microscopic reversibility holds. Property (iii) remains true by definition while (iv) continues 
to govern all interactions between the system and reservoir, or between different systems in an ensemble (Fig. [1}. 
We therefore find that the four conditions for detailed balance apply to a sheared fluid, amended by an additional 
conserved quantity, the total shear. These amended conditions give rise to a nonequilibrium counterpart to detailed 
balance, that can be derived from familiar statistical considerations in a straightforward way. 

We continue, below, to derive the exact consequences of the above conditions, in terms of the resulting microscopic 
transition rates. That derivation is performed without approximation, and therefore applies arbitrarily far from 
equilibrium. However, as with any investigation, the degree to which a real experimental system is approximated by 
the idealized assumptions depends on the fluid in question. For instance, the applicability of any exact equilibrium 
calculation depends on how well a real system has been equilibrated and how well its symmetries are represented by 
the theory. For the driven cases considered here, the question of how close the experiment comes to a true steady 
state is even less straightforward. In an experiment on an isolated system and reservoir, the fluid is in principle 
only in a quasi-steady state because energy is continually pumped into the system by the driving force, tending to 
heat the fluid. Nevertheless, the above assumptions can be experimentally realized to arbitrarily high accuracy. Let 
us elucidate the types of fluid to which the above assumptions apply. In many complex fluids a solvent acts as a 
thermostat for the system of interest, as it possesses many more degrees of freedom than the mesoscopic particles of 
interest, so the temperature remains relatively constant. In the limit of large ratio of solvent to complex degrees of 
freedom, the steady state condition becomes exactly realized. Since the relaxation time of the solvent is usually much 
smaller than that of the complex fluid, the solvent remains close to equilibrium even when the complex fluid is driven 
to a highly nonequilibrium state. In practice, experiments on complex fluids are able to obtain reproducible bulk 




(4) 



LOjiPj = WijJH, 



(5) 
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FIG. 1: A sketch of a large volume of sheared fluid: the fluid is hypothetically divided into a large number of subvolumes, 
the 'systems'. The height of each subvolume is much larger than the correlation length L corr such that interactions with the 
surroundings, i.e., the 'heat reservoir', only take place at the boundaries. The time evolution of a single system accumulates a 
certain amount of shear which is different for each system (indicated by the varying slopes). 



steady-state behavior, indicating that the system is insensitive to spatial or temporal temperature gradients. These 
form the class of systems for which NCDB is expected to hold, and that reproducible behavior is the subject of our 
exact theory. 

With these prerequisites in place, we can begin the derivation of NCDB from a nonequilibrium ensemble. Imagine 
a very large volume of fluid divided into a large number N of sub- volumes or 'systems'. These sub- volumes are still 
large enough that correlated regions within them are negligibly small (see Fig. [I]). The crucial idea in our derivation 
is to imagine this volume to be so large that it constitutes an ensemble in which the sub-volumes can be considered 
as representations of all possible realizations of our system under consideration. A certain constant shear rate J is 
applied to the volume by moving top and bottom boundaries (while perpendicular dimensions are infinite or periodic). 
Each system i then follows a particular trajectory in phase-space over a time period r accumulating a certain amount 
of shear. In order to ensure steady state properties we will eventually take r — > oo, so that any initial transient 
behavior has decayed. 

With the nonequilibrium ensemble constructed in this way, the probability distribution of trajectories can be found 
by Gibbs's familiar method for deriving the probability distribution of a large collection of countable objects, which 
is exact when no correlations exist between those objects, as is the case here. The probability that a system follows 
phase-space trajectory T is p(T) = nr/N, where nr is the number of times T is realized in the ensemble. Here, we take 
phase space and time to be discretized to make trajectories countable. Eventually the discrete intervals will vanish 
in the continuum limit and probabilities p(T) will be replaced by distributions p[T] T>T. 

The statistical weight of the ensemble Wn is the number of distinct ways we can arrange the systems in the 
ensemble, which is given by the usual combinatorial formula 

AM 

Wn = r, (6) 

Zjrnr' 

since systems following the same trajectory are indistinguishable. The most likely distribution of trajectories (the one 
adopted by the overwhelming majority of such ensembles) is the one with maximal statistical weight and can be found 
by maximizing the corresponding ensemble entropy Se — hi Wn- Applying Stirling's formula In AM w A" In N — N for 
large N yields 

Se = -N^2p(T)hxp(T). (7) 
r 

Thus, the ensemble entropy per system is the path entropy 



Sr = -$>( r ) ln P( r )- 
r 



(8) 
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This path entropy is familiar from approaches to nonequilibrium statistical mechanics in the spirit of Jaynes's method 
of maximum entropy inference (MaxEnt). In fact, the first derivation of NCDB was based on this information 
theoretic approach [3]. In the present framework of the nonequilibrium ensemble (cf. we are able to avoid certain 
controversial aspects of MaxEnt. Instead of a subjective interpretation of probabilities as representations of our state 
of knowledge, we have defined probabilities in terms of countable physical realizations. Furthermore, we have a clear 
idea of the prior set of trajectories that we reweight in the maximization procedure; it is the set of all physical 
trajectories with the same weights that they have at equilibrium. All unphysical ones are given zero weight a priori. 
The distribution p(r) is found by maximizing Sr subject to the constraint 

]>>(r)7(r) = Jr ( 9 ) 

r 

due to conservation of shear. Here, 7(T) is the total shear acquired by an individual system following path T. With 
J = 0, the maximization just returns p eq (T), the equilibrium distribution of trajectories defined as the prior for the 
calculation. In the driven steady state we obtain instead the result 

/ r (r) oc p eq {T)e^ {T) , (10) 

where v is the Lagrange multiplier associated with the shear constraint Eq. ^j. 

This form of the nonequilibrium distribution is a direct consequence of the way in which the system is driven, namely 
by fluid at the boundaries such that the equations of motion of the driven system are the same as in equilibrium. In 
fact, the measure defined on trajectory space for the driven ensemble is the same as in equilibrium, the trajectories 
are only reweighted under the additional shear constraint giving rise to the factor e vl ^ T \ 

Now that we are in possession of the statistical weight of nonequilibrium trajectories, we can consider the implica- 
tions of Eq. (|10[) on the individual transition rates, which give rise to all the complexity of the system's evolution. It 
is with this change of perspective that a useful and testable result can be obtained. 



III. RULES FOR TRANSITION RATES 



From now on we denote a transition rate between states i and j in the driven steady state with in order to 
distinguish it from the corresponding equilibrium transition rate uiij that satisfies equilibrium detailed balance. The 
transition rate is defined as the probability of making the transition i — > j per unit time At, where it is understood 
that At is so small that only one transition can occur. This means that transitions from i to j via a third state 
are neglected. In probability theory the definition of the transition probability P(j\i) is usually in terms of a joint 
probability P(j,i), i.e. P(j\i) = P(j,i)/P(i). If the distribution of trajectories p(T) is known, these probabilities of 
individual microstates are determined by counting those trajectories containing state i at time t and state j at time 
t + At. A transition rate is then defined as fly = P(j\i)/At in the limit of vanishing At, where 

p(M = E ^ i)er ^ } . (ii) 

The sum YlrlUj) denotes a summation over all trajectories T containing the transition i — > j. Using Eq. (|10[) a 
transition rate in the driven steady state can be written as 

Uij - At ~ Ai£ r|ier ^ ( r) e -Kir M 

Inserting unity as the delta-function J £(7 — 7(r))d7 yields 



SZo e" 7 Er|fe) g r % - l(r)) P e "(T) d 7 



We notice that the average over the delta-function 5(7 — 7(r)) is related to the conditional probability distribution 



(14) 
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The distribution p% q {"~l\i , j) contains the probability that the system at equilibrium accumulates a total amount of 
shear 7 over a time period r due to equilibrium fluctuations, given that it made a transition from i to j. Similarly, the 
denominator in Eq. (j 1 2 j) leads to an expression containing p% q {^f\i). The subscript r denotes the implicit dependence 
of these distributions on the duration of the trajectory T. The driven rate now reads 

Uij r l ™ AtE r|ier ^(r)r oo ^( 7 |z)e-d 7 • (15) 

Here, the r — > 00 limit guarantees that the system has attained its stationary state. After factoring out the equilibrium 
rate 

„ a . W>"f> (16) 

we can write the driven transition rate as 

X!° P% q {l\hi)e vl d 7 
Qij = uj tJ km gg . (17) 

We refer to this result as the canonical-flux representation of NCDB for a transition rate i — > j in the driven ensemble 
Note that the driven rate is a function of the Lagrange multiplier v, the flux conjugate parameter. This parameter 
is the analogue of the inverse temperature [3 which characterizes the heat transfer in canonical equilibrium systems. 
One notices that the driven transition rate is proportional to the equilibrium rate, enhanced or attenuated by a 
factor which is determined by the equilibrium statistics. If an equilibrium transition from i to j is likely to lead 
on to a certain amount of shear in the future, then this transition is enhanced in the driven steady state. As in 
the equilibrium case m/2 constraints act on m driven transition rates, therefore Eq. (|17[) can be considered as a 
nonequilibrium counterpart to detailed balance. The important conclusion is that for the class of driven steady states 
obeying the principles (i) — (iv) we obtain the microscopic dynamics of the nonequilibrium steady state from the 
corresponding equilibrium ensemble. 

Depending on the physical system, many set of equilibrium (prior) rates are possible, from exact and deter- 
ministic ones obtained from a Hamiltonian formulation to stochastic ones. Transitions which have zero probability in 
equilibrium (i.e., unphysical ones) are forbidden in the driven system as well. 



A. Alternative representations 

In the above derivation of NCDB the nonequilibrium ensemble was constrained by a fixed mean shear Jr (Eq. @). 
In analogy to equilibrium ensembles we can therefore consider Eq. (|17p as a canonical representation of the driven 
ensemble: shear is a quantity which can be exchanged between the system and the nonequilibrium reservoir with 
a fixed average. Indeed both numerator and denominator in the enhancement factor take the form of an average 
where the conditional probability p% q {"f\i, j) (respectively p% q {^\i)) is weighted with the shear distribution e" 1 over 
all possible shear 7. A microcanonical-flux ensemble allows only for a fixed total shear 70 . One can intuitively argue 
that in this ensemble the flux distribution has to be substituted by a Dirac delta function: e ul — > ^(7 — 70) (for an 
alternative derivation see Q). With this substitution Eq. (fTT|) reads 



fiy=Wy km eq . (18) 



This representation proves favorable when the conditional probability distributions of 70 are explicitly known. Similar 
to the ensembles in equilibrium statistical mechanics, microcanonical- and canonical-flux ensemble yield the same 
result in the thermodynamic limit, which is here the limit of large N and r. 

From Eq. (TTT|) it is also possible to derive a r-independent representation of the driven transition rates [4|. We 
reproduce the derivation here in full as it is used in the following investigation of NCDB in the context of discrete 
state spaces. Let us define the quantity Axji which denotes the shear contribution of the transition i —> j in time At. 
Using Axji we can write the conditional probability distribution p% q (^\i, j) as 

P?{-Y\i,3)=PrL At h-Ax ji \j), (19) 
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i.e., p^ 9 (7|i, j) is given by the probability distribution of accumulating the remaining shear 7 — Axji in the remaining 
time t — At, starting from state j. Defining the function iritis, r) as 



rrn{v,T)=hx p^( 7 |i)e^d 7 , (20) 

J - — OC 

allows us to rewrite Eq. (fT7|) in the form 

ln% = Um (in [°° p e T 9 _ At (j- Ax ji^e^dj-rmip, t)). (21) 

A change of the integration variable then leads to 

In— ^ = uAxji + lim [mAv, r - At) - mAv, r)] (22) 
= vAxji + lim [m;(i/,T) — mAv, r)l — C(^7 , (23) 

T— >00 

where 

C(^Ai) = lim [m,-(i/,r) - ra^{v,T - At)]. (24) 

T — >00 

An important property of the function £ is its state-independence, which follows from Eq. (|23p upon changing r — > 
r + At. To first order in At we then have 



d 

C{v,At)= lim —rrij{v,T)At. (25) 



We now introduce two important r-independent quantities [4]. Firstly, we identify 

Aq jt = qAv) - qi(v) = lim \mAu,T) - m t (v,T)}, (26) 

r — >oo 

and secondly, the state- independent rate of change of m,i{y, t) in the long-time limit is denoted as 

Q(y) = lim ^-m t ^, T ). (27) 

We therefore see that Eq. (|2"3")) gives rise to 

n tf (1/) = e ^,+A,,,(W-QHAi i (28) 

which is a r-independent representation of NCDB. Equation l|28p reduces to Eq. |(3J) in the continuous time limit, 
At^O. 

If we compare Eq. (|28p with Eq. (jTTJ) we realize that the ratio of equilibrium Green's functions is here translated 
into three distinct factors. The factor e uAxii measures the direct flux contribution of a transition and is large if either 
the flux carried by the transition i — ► j is large {Axji 3> 1) or the system is strongly driven (y 1). By itself 
this factor would simply boost every transition in the flux direction irrespective of the state space structure. The 
important extension to such mean-field ideas is expressed in the factor e Aqji ^\ where Aqji is formally defined as (cf. 
Eqs. (HU and P61) ) 



Aqji(v) = lim 



n /- OO oo^(7l0e^d 7 



(29) 



This quantity measures the increase (or decrease) in probability that the system will go on to exhibit the imposed 
shear 7 if it performs the transition i — ► j. Thus the rate of a transition not only depends on the immediate flux 
contribution, but also on the prospect for future flux. 

Formally, the function Q{v) can be defined as the scaled cumulant generating function of the shear current J = 7 /t 

Q(v)= lim - In / p e T q (J\i)e urJ dJ, (30) 

r^oo T J_ oc 
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FIG. 2: Example diagram of a basic graph for a five-state system. The dotted line denotes an exterior edge connecting the 
same type of states in different periods. Here states and 5, as well as states 1 and 6 are of the same type. 



which agrees with Eqs. PD|) and (f!?T)) upon transforming from the (extensive) shear 7 to the (intensive) shear current J. 
The distribution p% q (J\i) contains the probability to observe a shear current J over time r in the equilibrium system, 
given it was initially in state i. Eq. ([50)) expresses the fact that v and J are conjugate quantities in the sense of large 
deviation theory and implies that Q{v) is the Legendre transform of the rate function of J 0, Q. Consequently, Q(y) 
and J are related via 



In the following we refer to the function Q{v) as 'flux potential' due to the analogy with the usual thermodynamic 
potentials. 

In order to fully determine the transition rates of a system in the driven steady state according to Eq. ([28]) one 
therefore has to know the set of equilibrium transition rates, their individual shear contribution Ax™ as well as the 
flux potential Q and the functions Aqji. At equilibrium the are usually strongly constrained by detailed balance 
and symmetry considerations. The Axji are local and constant properties of the states. The crucial and non-trivial 
task is to find the potential Q and the set of Ag^-'s, which depend on the global structure of the state space. 



Master equation systems can be discussed intuitively with the help of a graph representation, in which vertices are 
assigned to the different states i of the system and edges to the possible transitions @. If a transition is physically 
allowed to take place, i.e., > 0, then equilibrium detailed balance demands that the reverse transition u>ji is also 
non-zero. Only connected graphs are considered in order to satisfy the requirement of ergodicity. For the following 
discussion of driven steady states further assumptions are necessary. In order to guarantee that the system can exhibit 
a constant macroscopic steady state shear current J, it is assumed that the state space has a periodic structure along 
a direction x, which is associated with the amount of shear accumulated. This assumption is made without loss of 
generality since one period may be arbitrarily large. A transition i — ► j contributes the shear increment A 
The periodic structure implies that, for a given period n, the j + nth vertex is of the same type as the jth with re- 
position shifted by a fixed amount. Likewise, if vertices i and j with 1 < i, j • < n are connected by an edge so are 
vertices n + i and n + j. 

We define the basic graph as the graph corresponding to the non-periodic connected set of n vertices. Its set of 
edges will be denoted interior edges in order to distinguish them from exterior edges connecting vertices of the basic 
graph with vertices of the next or preceding period. More precisely, an edge is exterior if it connects a vertex i of 
the basic graph with a vertex I where I > n + 1 or I < 0. Obviously, exterior edges always occur in pairs connecting 
the same states in different periods. Using this convention one has to bear in mind that two states of the same type 
can be connected by more than one edge, namely by one interior edge and multiple exterior edges (connecting, e.g., 
to the previous or next period, or both). See Fig. for a depiction of such a basic graph. The total number of 
transition rates in the system is J2?=i di, where di is the degree (or connectivity) of the ith vertex of the simple graph 
including exterior edges. The minimal number of transition rates in an n-state driven system is 2n which corresponds 
to a graph in the form of a simple connected path. For this class of state spaces the problem of finding the driven 
transition rates has a particularly straightforward solution (see Sec. IV Bl below). 




(31) 



IV. GRAPH REPRESENTATION FOR NCDB 
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FIG. 3: Example diagram of a three state system mapping onto a periodic graph structure. 



This notion of a periodic graph structure is basically a convenient way to visualize the current in the state space. 
Systems with a limited number of states would usually be depicted as a basic graph alone without any exterior edges. 
Yet, if the system exhibits a nonequilibrium steady state with some kind of current, the periodicity automatically 
arises as a consequence of the flux carrying transitions. Consider for example the three state system in Fig. [3l The 
only way the system can be in a driven steady state is by featuring a rotational current which is measured by the 
windings performed in time t. This then naturally maps onto a periodic network with a simple connected path as 
basic graph. The loop 1^2^3^1in the original depiction then becomes an external loop 1 — > 2 — > 3 — > 4, where 
state 4 is of the same type as 1 yet distinct due to the accumulated integrated current of one period. 

At this point it is appropriate to briefly discuss the distinction between external and internal loops, which is 
connected to Kolmogorov's criterion, an equivalent statement of equilibrium detailed balance [10 . 1 1 11 ] . In order to 
elucidate this we focus on transition rates in continuous time with At — > such that Eq. (|28[) assumes the simpler 
form of Eq. ©. 

If we consider a closed internal loop in the basic graph, i.e., a closed path leading back to the identical state such 
as, e.g., the loop 2-^3^4^2in Fig. we see that the product of transition rates following this loop in a given 
direction of rotation is the same as in equilibrium 



^23^34^42 — ^23^34^42 • 



(32) 



This is a simple consequence of the fact that both Axji and Aqji are given as differences of state properties, implying 
that along a closed internal loop the product of exponential factors vanishes. Therefore, for any internal closed loop 
in the basic graph we have 



^12 ' ' ' ^n-l,n^nl — ^12 ' ' ' W„-1„W„1. 



(33) 



Since the equilibrium rates satisfy detailed balance, the ratio of forward and backward transitions is given as W12/W21 = 
e -p[Ea—Ei) t it is then easy to see that the following equality holds 



Wl2 ■ ■ ■ U> n -i, n lLl n i 
^21 ' ' ' W n n _lCJl n 



= 1. 



Due to Eq. (|33p the same relation is true for the product of driven transition rates 
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• S2n-l,n"nl — ^21 ' ' ' fin,n-l"l 



(34) 



(35) 



Kolmogorov's criterion now states that equilibrium detailed balance holds if and only if Eq. (|35[) is satisfied for every 
closed path in state space. NCDB according to Eq. ([2]) satisfies this criterion for every internal loop. However, this 
does not lead to a contradiction, since Kolmogorovs criterion is violated for external loops in state space. For an 
external loop we have instead of Eq. ([35)) 



12 ■ 



^21 ' ' ■ ^n,n-l"l 



1' e 2vAx lfl 



(36) 



where V denotes the state of type 1 in the next period and Axyi is the accumulated shear or 'length' of a period. 
The conclusion of this discussion is that the graph representation for NCDB in terms of a periodic graph structure is 
consistent if the system does exhibit a steady state current. In this case NCDB leads to a violation of Kolmogorov's 
criterion as expected. 
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At equilibrium the validity of detailed balance, or Kolmogorov's criterion respectively, implies that for every two 
states k and I the ratio 

_ U>kiU>i2 ■ ■ ■ u> n ^i in u) n i 
■Jkl = (o< J 

is independent of the path between k and I. Out of equilibrium this path-independence is not generally expected. 
Instead one can consider InSu as an 'action functional' associated with a particular path and derive a fluctuation 
relation for the entropy production [l3j . 

V. TOTAL EXIT RATE RELATION 

In the remainder of this chapter we focus on NCDB in the context of continuous time Markov chains, where NCDB 
is expressed in the form of Eq. @. In this case it is possible to derive a fundamental relationship between the flux 
potential Q{v) and the driven transition rates which is a central result of this article and leads to a variety of important 
conclusions for NCDB. 

Consider an individual state i connected to di other states. The system in state i spends a random time until it 
makes a transition to one of the connected states. For a continuous time Markov chain this waiting time has the 



exponential distribution [ill ] : 



h i (t)=* i e-"* t , (38) 

where the total exit rate is defined as cr^ = The probability that the particle jumps to site j is then 

Pij = LOijj <Ji. Our quantities of interest are the conditional probabilities or Green's functions p% q {^f\i) which, if known, 
would fully specify the driven transition rates via Eqs. (|29p and The Green's function for state i can be determined 
by the following considerations. From state i the system can only perform a transition to a connected state j within 
the network, from where its further displacement is determined by the Green's function of that state j. Taking into 
account the waiting time in state i and the probability to perform the jump to state j, p% q {~f\i) is therefore related 
to the di Green's functions of the neighbouring sites according to 

p?( 7 |i)= f T dth^r-t^P^PtHl-^^ + M^Sh), (39) 

where tpi (r)dr denotes the probability that no jump has occurred 1 out of state i up to time r: ipi (r) = 1 — hi (t)dt = 
g-wiT j n £ ne next step we introduce the quantities rrii(v, r), Eq. (|20[) , which are ultimately related to the functions 
qi{v) in the long time limit via Eq. (|26[) . Multiplying Eq. (|39[) by e vl and summing over all possible shear 7 from 
state i yields, after a shift in the summation variable 7, 

e m, ("- T ! = f di/i,(T-f)VP !3 e m ^^ +1/Al "+e-''' r . (40) 

In the next step we can substitute Pij, and hi(t) (Eq. (|38|) ). Rearranging terms then yields 

e m '("' r ' +fffr = f dte aa Vw i je mj( "' t)+l ' Alii + 1. (41) 

The integral can be removed by taking a derivative with respect to r on both sides. This leads to 

^mi(i/, t) + <r^J e ™d^T)+* tT = e a lT uj i:j (Mlvri+vAxn ( 42 ) 

Or, after further rearrangement, 

A m .(^ r ) +(T . = y U0 . e mj(",T)-m,(i',T)+»Ax J(i (43) 

h} 



1 The contribution of ipi(r)5(-f) to the Green's function p% q ("j\i) was erroneously neglected in the derivation presented in Ref. [fj. Never- 
theless, the final result Eq. 1451 remains unchanged. 
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In the long time limit Eqs. ((26)) and (l27l) hold. The result of taking the r — > oo limit in Eq. (|43f is thus 

Q(i/) + o-i = ^M+i'A-ii. (44) 

3 

On the right-hand side we can identify the transition rates in the driven steady state according to Eq. |(3J). We 
therefore obtain a fundamental relationship between the equilibrium transition rates, the corresponding rates in the 
driven steady state, and the flux potential Q(v): 

QO) = -a;, (45) 

where Ej(i/) = Eq. ([43]) states that, for every state i, the total exit rate in the driven steady state differs 

from its equilibrium counterpart only by a flux dependent (but state independent) constant. On the basis of this 
central result a number of important implications of NCDB can be derived. It turns out that it is not necessary to 
know the Green's functions of the equilibrium system in order to determine the Ag's and the driven transition rates. 
Rather, as we will see more explicitly below, the quantities of the NCDB formalism are intrinsically related to the 
graph structure via Eq. (|4"S"|) . 

A. Invariant quantities 

It is now straightforward to formulate two sets of invariant quantities for the sheared steady state. The first was 
found in Ref. |4| and is a consequence of the asymmetric property of Ag^ and Axji. From Eq. the 'product 
constraint' follows directly, 

QijQji — (46) 

Secondly, Eq. (|45[) directly implies the 'total exit rate constraint', 

- Ej = o-j - o-j. (47) 

We therefore find that in the driven steady state both the product of forward and reverse transition rates and the 
difference of total exit rates for every pair of microstates are the same as in equilibrium and therefore invariant with 
respect to the driving. No near-equilibrium assumptions have been made in the derivation, so the above relations are 
both exact and valid arbitrarily far from equilibrium. 

With the formalism devised here, the task of finding the microscopic nonequilibrium dynamics given an imposed 
macroscopic current is greatly simplified and follows straightforward rules. The above relations are furthermore 
accessible to verification in an experiment or simulation [l2[ • 

B. Network rules 

With these results we are able to devise a systematic method to determine the driven transition rates for arbitrary 
graph configurations. Considering the whole basic graph with n states, there are n equations in the form of Eq. 1(15]). 
This set of equations is sufficient to determine all the unknown quantities, i.e., the flux potential Q and the Ag's, as we 
now show. Due to the relationship Ag^ = — Aqji, every edge of the basic graph is associated to two transition rates 
(forward and backward transitions) depending on one Aq. The number of independent Ag's is further constrained by 
closed paths (loops) in the graph, i.e., paths that begin and end at the same type of vertex, because the sum of Ag's 
along such a path is obviously zero (loop constraint) . One can then easily see that the total number of independent 
Aq's in the basic graph is always n — 1, as follows. Consider first the most simple basic graph configuration, namely 
all n states connected as a simple path without any loops. In this case the number of edges is trivially n. Since one 
loop constraint is generated by the periodicity, there are n — 1 independent Aq's. From this simple connected graph 
all graphs of higher degrees can be generated by adding new edges. But adding an edge generates a new Aq and at 
the same time a new loop constraint, so that the number of independent Aq's always remains n — 1. 

We formulate the following network rules for the calculation of the driven transition rates in networks of arbitrary 
connectivity: 

• Edge rule. Every interior edge and every pair of exterior edges in the basic graph corresponds to two rates 
containing the dependence on one Ag. The driven transition rates are given by Eq. 
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• Vertex rule. For every vertex in the basic graph the difference between the driven and equilibrium total exit 
rates equals the flux potential Q (Eq. P5jl ). 

• Loop rule. For every closed path of edges the sum of the Ag's along this path is zero. 

In this formulation there are in total n equations and n unknowns, namely one Q and n — 1 Ag's. The number of 
independent equations can always be further reduced by eliminating Q, such that one is essentially left with n — 1 
equations for n — 1 unknown Ag's. The solution of this system of equations fully specifies all the driven rates in the 
system as well as the flux J which is related to the flux potential via dQ(v)/dis = J, Eq. (f3Tj) . In the given framework 
the driven transition rates and the flux are determined as functions of the flux conjugated parameter v. 

An alternative, ^-independent representation is based on Eqs. (|46j) and (|47| . However, the number of constraints is 
here not sufficient to determine all the driven rates for general network structures. On the one hand there are X)T=i 
product constraints and n — 1 exit rate constraints. On the other hand, for an arbitrary graph configuration, there 
are Y^i=i ^ transition rates. Therefore only for graphs with the topology of a simple connected path (where di = 2), 
we can determine the rates completely from the invariant quantities without using the network rules. In this case 
we have 2n transition rates and 2n — 1 constraints stemming from the exact relations. The transition rates are 
fully determined if additionally the relationship between the transition rates and the current is provided. In this 
formulation the driven rates depend on J directly instead of being parametrized in terms of the parameter v. The 
relationship between current and rates is further elucidated in the next section. 



VI. CURRENT RELATIONS 



In Master equation systems the stationary current J is generally defined as the average of the flux contributions of 
every state i over the steady state distributions pi 



J 



{./} 



At equilibrium J is identically zero due to detailed balance, Eq. ([5]), and the property Aa 



(48) 



-Axij. In general, 



at steady state, the distribution pi is the solution of the Master equation (J3J) under the condition of station ary 
dpi(t)/dt — and normalization YliPi = 1- The solution can be formally obtained by 
In matrix notation Eq. ([3]) is given as (writing the transition rates now capitalized): 





-Si 


n 2i 


n 3 i • 




\ 


/Pl\ 






-s 2 


^32 ■ 


• Sln2 




P2 






^23 


-S 3 • 


' ^n3 




P3 


V 




^2n 


^3n - 




/ 


V Pn I 



0, 



(49) 



i.e., as MjPj = 0, where M is a n x n matrix with entries 

My = % - fyEi. 



(50) 



It has been noted earlier that, due to the periodic graph structure, in addition to the internal edge, there can be 
multiple exterior edges connecting the same types of states in different periods. These additional rates are contained 
in flij which denotes the sum of all transition rates from a state of type i into a state of type j. In the Master equation 
this is considered in the summation {j} over the set of adjacent sites which includes the exterior edges. Denoting the 
transition rate matrix including the normalization (e.g., in the first row) by M, that is 



M 



/ - 1 1 - 1 

f2l2 — £2 ^32 
^13 ^23 



V f^ln ^2 



3/2 



- 1 \ 

^n3 

-£„ J 



(51) 
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the Master equation including normalization reads 



/ - 1 



1 



_ 12 

fil3 ^23 



V Qln ^2n fisn 



1 • 


• 1 


\ 


f P l\ 




fi 3 2 • 




P2 


(°) 


-E 3 • 


■ fin3 




P3 


= 


^3n ' 


' — E n 


/ 


\pn) 


W 



(52) 



We therefore find the occupancies pi by matrix inversion as M _1 ei, i.e. in the first column of the inverse matrix M _1 . 

A new relationship between the current J and the transition rates is obtained when we take the derivative (denoted 
by a prime) of Eq. (|45| with respect to v. 



g' = ^(A^ + A^)^, 
{j} 

for every state i. Here and in the following the ^-dependence is dropped. Using Aq'^ = q'j 
equations as 

Q' + q'jT.i - qfaij = Axjiflij. 



(53) 

q[ we can rewrite these 
(54) 



Since all relevant physical information is contained in the differences rather than the individual q^s, we have the 
freedom to fix one boundary value. If we choose q\ = const, obviously q[ — and thus q[ = —Aq'^. Defining the two 
column vectors q' = (Q',Aq' 12 , Ag lT[ ) T and a = (ai, ...,a n ) T with a, = ^- x ji^iji we can cas t the system of 

equations (|54[) into the following matrix form 



(\ fi 12 Ol 3 
1 — S2 ^23 
1 O32 —S3 



V 1 ^. 



CI 



II?} 



filn \ 

fi2n 

fi3n 



/ \ 

V Agi„ y 



0.3 
V a n / 



(55) 



We realize that the left hand side reveals the transpose matrix of M. Eq. (|55|) is therefore equally expressed as 

M T q' = a. 



(56) 



With this result we can rederive the fundamental relationship Eq. (|3"Tj) between the flux potential Q and the flux J. 
Eq. (|56p is formally solved via matrix inversion 



l' = (M r )" 1 



a. 



(57) 



For any matrix the transpose of the inverse matrix is the inverse of the transpose matrix: (M T ) 1 = (M 1 ) T . Since 
M _1 ei is the formal solution of the Master equation, we can deduce that the first row of (M T ) _1 contains the steady 
state occupancies pi. Solving Eq. (|57| for Q' with the given expressions therefore leads to 



Q' = Xlp»y^ AxjjQ 
{j} 



(58) 



so that by comparison with Eq. l[48|) . we can conclude that Q' = J, i.e., Eq. (|31| holds. We have therefore shown that 

(59) 



J = J £(Aq , ji + Ax ji )Q 
{3} 
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for every state i, which is a new representation of the steady state current J in terms of Aq'^ instead of the microstate 
distributions pi . 

At this point it is instructive to associate a physical interpretation with the v derivative of the Ac^. Taking the 
derivative of Eq. ([29]) with respect to v leads to the following expression: 



Aq' 3l = lim 



Hg lP% q {l\3) ^ d 7 fZolPT q (.l\i) ^ d 7 



(60) 
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In the discussion of the path-entropy maximization in Sec. Mil we obtained the fundamental result (Eq. (JTOj) ) that in 
the driven ensemble the trajectories are reweighted with the exponential factor e vl . This implies that in the present 
context p% q (^\j) e" 7 can be related to the unnormalized conditional probability that the system exhibits shear 7 from 
state i under the driven dynamics, i.e., 

pf ( 7 H) oc p?fr\()er>, (61) 

in the limit of large r. Introducing 



no7P^(7|t)e^d 7 
L°° 0O Pr 9 (7K)e^d 7 



( 7i (r,i/)) = -jH — J q — — J _ , (62) 



which denotes the mean shear that the system accumulates over time r in the steady state from state i, one can then 
express Eq. (f6T)|) as 

A<^ = lim [(7j(r,v))-(li(T,v)}], (63) 

T — >00 

i.e., as the difference in mean shear between states i and j over infinite time. Furthermore, from Eq. (|57[) we find that 
the quantities Aq'^ can be determined using Cramer's rule [lOj : 



det 


"(M T )W 


det 







= 1 r _ 1 J , (64) 



where (M T )W means that the ith column of M T has to be replaced by a. From the set of Aq' u all other quantities 
Aq'^ follow by subtraction: Aq'^ = Aq' u — Aq[-. 

Relationship (|59p is remarkable because the stationary shear current is determined from the properties of a single 
state and its neighbors only, instead of the average Eq. (|4"5| over the whole basic graph. It does not involve energetics, 
since the knowledge of two constant 'shear values' in addition to the transition rates is sufficient to determine the 
current. One constant, Axji, measures the immediate difference in shear between states i and j and the other, Aq'^, 
the difference in mean shear that the system accumulates in the steady state over infinite time from states i and j. 
The set of Aq'^ is directly related to the transition rate matrix M via Eq. (|64[) . 

From a computational point of view, no advantage is gained by determining the flux via Eq. (|59p if the quantities 



Aq'^ are determined by the same formal matrix inversion method as the probability distributions. However, Eq. 
provides us with a new interpretation of the steady state current and of the elements of the inverse transition rate 
matrix M _1 . The quantities Aq'^ have a precise physical meaning independent of the NCDB formalism. It might be 
possible to identify related current expressions for other steady state systems that are not contained in the class of 
systems for which NCDB is valid. 

A. Fluctuation relation for the shear current of a trajectory 

Fluctuation theorems are mathematical relations for the fluctuations of thermodynamic quantities like heat (entropy 
production) or work in nonequilibrium systems (see [l4| and references therein). In a nonequilibrium steady state 
the so called asymptotic or steady state fluctuation theorem (SSFT) states that the probability distribution p T (p) of 
finding a particular value of the thermodynamic quantity p over time r satisfies a certain symmetry relation of the 
form (cf. El EE EE El) 

Pr(P) ^„cpr^ ( g 5) 



Pr(-p) 

where = indicates the asymptotic behavior for large r and c is a constant. The SSFT Eq. (|65p (also refered to 
as Gallavotti-Cohen fluctuation theorem) represents a refinement of the second law of thermodynamics in that it 
quantifies the probability of observing temporary second law violations (negative p) in the nonequilibrium steady state. 
Relations similar to Eq. (|65p have been derived for a variety of systems with different thermostatting mechanisms. For 
deterministic systems in a compact phase space the SSFT for the entropy production is expected to hold universally 
under the chaotic hypothesis [161 ]. However, stochastic systems do not exhibit the same generality in their fluctuation 
behavior. Here, the validity of Eq. (|65p relies, e.g., on the characteristics of the noise, or on the thermodynamic 
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quantity considered. In a particular paradigmatic nonequilibrium particle model the SSFT has been shown to hold 
for the heat fluctuations, when the system is thermostatted by an equilibrium heat bath with Gaussian white noise 
characteristics 18]. However, when one considers the heat fluctuations, or the work fluctuations under the influence 
of non-Gaussian noise, the SSFT is violated [H H3, HI [H ||| . 

In a shear flow, a fluid region in the bulk of the volume receives noise from the sheared fluid surrounding it, which 
represents a nonequilibrium heat bath. NCDB quantifies the stochastic influence of such a nonequilibrium heat bath 
on the dynamics of the fluid. In systems satisfying NCDB a relation in the form of the SSFT holds for the fluctuations 
of the shear current of an individual trajectory and can be derived in a straightforward way. Our starting point is 
the nonequilibrium path distribution p dr (T) of Eq. (fTO]) . Using this path distribution, the distribution of the shear 
current j of trajectories of duration t can be formally expressed as 

T,rP ( r ) 

where the total shear of a phase space trajectory T is denoted by 7(r). Substituting Eq. (TIT)]) yields 

Here, the numerator is just the unnormalizcd probability distribution of observing the shear current j at equilibrium. 
Clearly, this distribution has to be symmetric under a change of sign of j, i.e., at equilibrium the probability of 
observing a shear current j over time t is the same as observing —j over the same time period. Taking the ratio 
Pr(j)/pT{—j) then immediately implies 

Pr(j) „ 2»j T 



e lvJ \ (68) 



for large r. The distribution p T (j) depends implicitly on the parameter v. 

The current fluctuation relation Eq. (|68p characterizes the fluctuations of j in the steady state and implies that 
the system's trajectory is overwhelmingly likely to follow the direction prescribed by the imposed driving, which is 
specified by the sign of the flux conjugated parameter v. Similar fluctuation relations for the current have previously 
been derived in the context of lattice gas models [24[ (see also [l4| and references therein). In the present case we have 
shown the validity of such a current fluctuation relation for a very large class of nonequilibrium steady state systems, 
namely fluid flows under continuous shear, as described by NCDB. Eq. (f6"8")) follows then very naturally due to the 
Gibbs property of the path distribution p dr (T). The general relationship between the SSFT and Gibbs distributions 
has been discussed in [25| . 

In the remainder of this article we apply the formal results of the previous sections to two simple Markovian jump 
models. 



VII. ANALYTIC SOLUTION OF THE TWO STATE GRAPH 



As our first example we discuss a basic graph with two different types of states 1 and 2. For generality, an asymmetric 
separation between the states is assumed. Two distinct arrangements are then possible: (i) a simple connected path 
leading to a zig-zag shaped graph, and (ii) a graph in the form of connected triangles. (Other connectivities are also 
possible, one of which was considered in Q.) For both types of graphs the driven transition rates can be calculated 
exactly using the results of the previous sections. Obviously, (i) is a special case of (ii) where the horizontal transition 
rates are set to zero. For clarity we first discuss the simpler zig-zag graph in some detail in Sees. IVILAl and [VII Bl and 
then present results for case (ii) in Sec. IVII CI 

Let us first consider a zig-zag arrangment of the two states leading to a ratchet shaped state space. Although the 
solution of this two state model follows quite intuitively, we apply the network rules of Sec. IVBl in a systematic way. 



A. Network rules 



Edge rule. The basic graph in Fig. [6] reveals one interior edge and one pair of exterior edges. In the following we 
use an intuitive notation for (u)pwards and (d)ownwards rates in positive (+) and negative direction (-). Driven rates 
are capitalized while their equilibrium counterparts are lower case. To the interior edge the transition rates U + = fl\2 
and D~ = Q21 are associated, which depend on Aq2i and carrying the shear iAa^i- To the pair of exterior edges 
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FIG. 4: Basic graph of the two state zig-zag model. The dotted lines denote exterior edges. 



we associate the rates U~ = 17 10 and D + = tt 23 , which depend on Ag 32 and ±Ax32 (note that Aq 32 = Aq 10 and 
AX32 = Ax±o). At equilibrium the number of different transition rates remains four due to the asymmetry of the 
state space, but equilibrium detailed balance requires that u~ /d + = u + /d~ = e~' 3AB , with AE being the energy 
difference between states of type 1 and 2. Throughout the following calculations the inverse temperature (3 is set to 
unity. According to Eq. |2]) the driven transition rates are then given by: 

V - = u - e -vAx 3 2-Aq 3 ^ D + = d + e vAx 3 2+Aq 32 _ ( 69 ) 

Vertex rule. There are two distinct vertices in the basic graph. To each corresponds an equation in the form of 
Eq. gU) 



Q = U + + U~ - (u + +u~), 

Q = D+ + D~ - (d + +d~). (70) 

Loop rule. There is one external closed loop 1 — > 2 — > 3 corresponding to subsequent transitions U + and D + . The 
loop constraint for the Aq's then reads 

Aq 21 + Aq 32 = 0. (71) 

In total we have the system of equations 

u+e „A X21 +A q21 +u - e -vA X32 -Aq 32 _ ( u + + u ~) = ( 72 ) 
d - e -vAx 21 -Aq 21 + d+e uA X32 +Aq 32 _ ( d + + d -j = q ^ 

Aq 21 +Aq 32 = 0, (74) 

for the three unknowns Q, Aq 2 i, and Aq 32 . The solution can be found in a straightforward way. We simplify notation 
by setting Ag = Aq 2 i = —Aq 32 , as well as Airi = Ax2i and Ax 2 = Ax 32 . Furthermore, we use the notation for the 
total exit rates o\ = u + + u~ and a 2 = d + + d~ . Substitution of the loop constraint yields 

{u + e vAxi +u-e-' jAx2 )e Aq -en = Q, (75) 
(d-e- uAxi + d + e" Ax2 )e~ Aq - a 2 = Q. (76) 

Elimination of e Aq in this set of equations leads to a quadratic equation for Q, namely 

Q 2 + (tn + a 2 )Q + (axa 2 - (u+e^ 1 + u~ e^ Ax2 )(d- e^ Axi + d + e» Ax2 )) = 0. 

(77) 

Due to equilibrium detailed balance we furthermore have 

o-to-2 - {u+e uAxi + u-e- uAx2 ){d- e - uAxi + d+e uAx2 ) = 2tt+d+(l - cosh(^(Aa; 1 + Ax 2 ))). (78) 



17 



The positive root then reads 



Q(v) = -V(oi + a 2 ) 2 + 8u+d+(cosh(i/(Aa:i + Ax 2 )) - 1) - -(01 + 05). (79) 



(80) 



In turn, Ag is determined from Eq. (|75p as 

Q + CTl 



Aq = In 



For a given set of equilibrium rates, the driven transition rates are now completely determined as functions of the 
flux conjugated parameter v. With Q and Aq we can furthermore calculate the shear current J(f) = dQ(v)/dv and 
the mean shear difference Aq' — dAq{v)/dv (see also below). Equilibrium is characterized by v = 0. As expected 
the flux then vanishes: J(v = 0) = 0. One can also easily see that for v — both Q(0) — and Aq(0) = 0, and the 
driven rates Eqs. (|69[) reduce to the equilibrium rates. 

B. Invariant quantities 

Since the asymmetric two state zig-zag graph is obviously a simple connected path, we can alternatively determine 
the driven transition rates using the z/-independent representation. To this end we have to set up the invariant 
quantities and find the relation between the rates and the particle flux. The invariant quantities follow immediately 
from Eqs. (gBJ) and (j4"Tj) : 

U + D~ = u + d~, (81) 
D + U~ = d + u-, (82) 
D+ + D- - (U + + U~) = d+ + d- - (u+ +u~). (83) 

Finally, to find the flux, we calculate the probability distributions of states 1 and 2. The Master equation for the 
two state graph has the following matrix form 



-(U+ + U-) D+ + D- \ f pj, 
U+ + U~ -{D+ + D-) JU 



= 0. (84) 



Including the normalization in the transition matrix leads to 

, S ) = ( 

where 



M(!M=( n ), (85) 



Matrix inversion then yields the probability distributions (or occupancies) pi , p\ , which, in turn allow for the calcu- 
lation of the current J via Eq. (|48p . The results for the probability distributions as functions of the driven transition 
rates are 

D+ + D- 

Pl = U+ + U-+D++D- * (8?) 

P2 = FTFTF+F' (88) 

According to the results of Sec. IVII we can alternatively determine the current and the mean shear differences Aq' 
from Eq. (|57[) . This equation is given by 

J \ _ (^-i\ T ( U+A Xl - U-Ax 2 
-Aq' 
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FIG. 5: Basic graph of the two state loop model. The dotted lines denote exterior edges. 



The matrix algebra is straightforward and yields 

U + D + - U~D~ 

J = (A Xl + Ax 2 ) — =— ; , (90) 

y 1 2, U+ + U-+D++D-' V 1 

= (D+ + U-)A X2 -(U+ + D-)A X1 

U+ + U~ + D+ + D- ' v ' 

The last two expressions agree with the ones obtained in the ^-representation, if the driven rates are substituted 
according to Eq. ([69|l with the known expressions for Q and Aq, Eqs. (|79|) and ([80l . 

For a given set of equilibrium rates and a prescribed particle flux, the transition rates in the driven steady state are 
determined as solutions of the set of equations (fSTj) — ([55]) and (j9"0"|) . In this case the driven rates are parametrized by 
J. 



C. Two state graph with a loop 

If we add two additional external edges 2 — > and 2 — > 4 to the two state zig-zag graph of Fig. [6l we obtain a state 
space with the structure of Fig. [5J The system in state 2 can now choose to reach state 4 directly or by going first 
'downhill' to a state of type 1 and then 'uphill'. The choice will depend on which path is more favourable to achieve 
the imposed flux J for given parameter values. As before, NCDB precisely quantifies the change in the transition 
rates under driving. 

We denote the transition rates associated to the additional exterior edges with H + = fl 2 4 and H~ = f^o- Both 
rates are independent of Aq since they connect states of the same type and symmetry requires that the corresponding 
equilibrium rates are equal. The rates are therefore fully specified as 

H+ = he u{Axi+Ax2 \ H+ = he~ u{Axi+Ax2) . (92) 

The effect of the additional edges on the other rates can be calculated analogously to the zig-zag graph. The total 
exit rate relation for state 2 (used in the vertex rule) has to be extended by the new rates, whereas the relation for 
state 1 remains unchanged. This leads to 

Q = U + + U- -(u + + u~), (93) 
Q = D + + D~ + H+ + H~ - {d+ + d~ + 2h). (94) 

Introducing the total exit rate a 2 = d + + d~ + 2h and substituting the expressions for the driven rates, we can write 
the system of equations as 

(d- e - ,/Axi +d + e ,/Ax2 )e- Aq + 2hcosh(v(Ax 1 +Ax 2 )-cr 2 = Q, (95) 

{u + e vAxi +u-e- vAx2 )e Aq -o 1 = Q. (96) 

The quadratic equation for the flux potential then reads 

= Q 2 + (<7i + 02 — 2/icosh(^(Axi + Ax 2 ))Q + o\o 2 — a\2h cosh(^(Airi + Ax 2 ) 

~(u+e» Axi + u-e- uAx2 )(d-e- vAxi + d+e» Ax2 ), (97) 
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(a) J (b) 

FIG. 6: (a) The four transition rates of the two-state zig-zag model (h=0.0) plotted as functions of the current J. (b) The 
function Aq(v) of Eq. (|8U[I (with Q(v) given by Eq. (|98[l 1 for two values of the equilibrium rate h. Parameter values: AE = 2.0, 
d+ = 1.0, d~ = 0.8, Azi = 1.0, Ax 2 = 0.5. 



with the positive root 

Q(v) = - y/Jal +a 2 - 2/i cosh(i/(Axi + Ax 2 ))) 2 + 8{aji + u+d+)(cosh(i'(Ax 1 + Ax 2 )) - 1) 

~2^ ai + ° 2 ~ 2h cosh (K Aa; i + Aa; 2)))- (98) 

As above, Aq is determined from Eq. ((80|) . 

Having thus specified all the driven rates, we can determine the remaining quantities in the matrix formalism. Since 
the rates H + and H~ connect two states of the same type, their effect on the occupancies balances to zero. The 
transition rate matrix M is therefore identical to Eq. and the probability distributions are given by Eqs. ([57]) 
and ([55)1 as in the zig-zag case. The shear current J and the mean shear difference Aq' are determined by Eq. (|57[) . 
which reads here 



( 3 \-(ik- x \ T { U+Axx-U-Ax 2 
\-Aq')-\ m ) \{D++H+ -H-)Ax 2 -{D- +H~ - H+)A Xl 



(99) 



We thus obtain 



j </k , U + D + - U-D- + (U + + U'){H+ - H~) 

J = (Ax 1+ Ax 2 ) D++D-+U+ + U- ' (1 ° 0) 

, = (D+ + U-)Ax 2 - (U+ + D-)A Xl + (H+ - H~)(A Xl + Ax 2 ) 

q D+ + D-+U+ + U- ' 1 ' 

In the limit h — > the two-state loop model obviously reduces to the zig-zag model discussed in the previous section. 
The four transition rates U + , U~, D + , D~ in this case are plotted in Fig. HJa) parametrically as functions of the 
current J for given values of the parameters AE, d + , d~ , Ax%, and Ax 2 . Clearly, the transition rates in the direction 
of J arc enhanced for larger J, while the transition rates in the opposite direction are attenuated. For large J the 
transition rates in the direction of the driving become proportional to J, as implied by the current relation Eq. (|100[) . 
The occupancies p\ 1 p 2 , which at equilibrium (J = 0) are determined by Boltzmann's law, converge to the value 1/2 
in the limit of J — > oo, indicating that the energy difference between states 1 and 2 becomes irrelevant for strong 
driving (see Fig. 0a)). 

When h ^ the system is able to gain a shear increment Ax\ + Ax 2 by going directly from state 2 to state 4, i.e., 
another state of type 2, without first going to state 3, a state of type 1. This implies that, under (strong) forward 
driving, it is advantageous for the system in state 2 to make the transition 2 — > 4 instead of the transition 2 — > 3, 
because it is then able to gain an increment Axi + Ax 2 in one step instead of two. This intuitively expected behavior 
of the two-state loop model is evident in Fig. [HJ where we plot the three forward transition rates under forward driving. 
While in the zig-zag case the rate D + is enhanced under increased forward driving, it is attenuated in the loop case 
due to the presence of the additional edge. The system finds that being in state 1 is less favorable in order to achieve 
shear, than being in state 2. 

The quantity Aq = Aq 21 , Eq. (SO]) (with Q of Eq. ([5S])). is plotted in Fig. E^b) for the two cases h = and h ^ 0. 
In the zig-zag case (h — 0) one notices that, for the given choice of the parameters Axi > Ax 2 , Aq becomes negative 
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FIG. 7: The occupancies pi and p2 of Eqs. (|87[1 and (|88[1 plotted as functions of the current J for parameter values AE = 2.0, 
d + = 1.0, d~ = 0.8, Axi = 1.0, A12 = 0.5. (a) In the two-state zig-zag model both distributions converge to 1/2 for large J. 
(b) In the two-state loop model, with h — 0.5. The occupancy of type 2 states, pi converges to 1 for large J, while pi decays 
to zero. 




for large forward driving (large is), indicating that state 1 has a larger propensity for future shear than state 2. This 
is because from state 1 the system can gain a larger (forward) shear increment (Axi) than from state 2 (Aa^)- F Q1 ' 
weak driving (small v) a crossover behavior is observed, where Aq exhibits local maxima and minima. In the loop 
case (h ^ 0) state 2 is favored due to the additional edge and thus Aq increases for both larger forward and backward 
driving. In fact, from the analytical expression Eq. (|80|) one finds that Aq(v) oc \v\ for large ±v in the case h ^ 0, 
while Aq(v) oc v(Ax2 — Axi) in the case h — 0. 



VIII. THREE STATE GRAPH 



A basic graph with three distinct states allows for a variety of periodic graph structures with different connectivities. 
In the following we consider a basic graph that contains an internal closed loop (see Fig. [9]). For this arrangement the 
driven transition rates can only be found using the network rules, since the number of independent rates exceeds the 
number of invariant quantities plus current relation. 

Edge rule. There are three internal edges and one pair of external edges in the basic graph. The three internal 
edges correspond to six distinct rates: f2i2, ^13, ^23, and the associated reverse rates. The pair of external edges is 
associated with i7io and ^34, where state is of the same type as 3 and 4 of the same type as 1. The set of driven 
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FIG. 9: Basic graph of a three state model with an interior loop. The dotted lines denote exterior edges. 



transition rates is specified as in Eq. J2]| 



^ij ^ij e J ^ 



Vertex rule (Eq. I|45p) gives one relation for each vertex. These three equations read 

Q = O12 + fil3 + Oio - (u>12 + U>13 + U>w), 

Q = ^23 + ^21 - (^23 +^2l), 

Q = n 31 + n 32 + ^34 - (^31 + W 3 2 + W34). 

3 and one internal loop 1 — > 2 — > 3 



Loop rule. There is one external loop 
Ag's read 



1 



Aq 2 



(102) 
(103) 
(104) 

1 so constraints on the four 



Agio 
Aq 32 



Aftsi = 0, 
Agi 3 = 0. 



Furthermore the internal loop requires that Aa; 31 = Ai 2 i + Aa; 32 . We use the following simplifying notation: 



In z 


= Aq 2 i, 


lny 


= Aq 32 , 







The loop constraints then imply that Agoi = Aq 3 i = ln(yz). In this notation the vertex rules read 



(W 31 

where the total exit rates are as before Ui 



W12 z + (W10 + Wi 3 )yz -0-1 = Q 

W23 y + W21 z~ x 
-1 



0-2 
0-3 



W 3i )(yzr 1 + W32y- 1 
Y], LOij . Elimination of Q then leads to two quadratic equations 



(W w + W 1 3)z 2 y + W 12 z 2 -W23yz + (a 2 -a 1 )z-W2i = 0, 
W 23 y 2 z + (cr 3 - a 2 )yz + W 2 i y - W 32 z - (W31 + VF 34 ) = 0. 



(105) 
(106) 



(107) 
(108) 
(109) 



(110) 
(111) 
(112) 



(113) 
(114) 



One realizes that already for the three state model an analytical solution is exceedingly difficult to obtain. Solving for 
example Eq. Q113P for y and substituting into Eq. (|114[) yields an equation for z which is of fifth order and therefore 
not exactly solvable. Alternatively one can determine the driven transition rates in a straightforward way by solving 
Eqs. (|1 13[) and (|1 14[) numerically for y and z. 

The current and probability distributions are easiest obtained from the matrix methods outlined in Sec. IVII For 
the three state graph of Fig. [5] the master equation reads in matrix form 



-(ft 13 + ^i 2 + 
^13 + ^10 



n 



10 ) 



^21 

-(Q 2S + n 2 i) 

«23 



0,31 + Q34 
^32 

-(fi 3 l + ^32 + ^34 




(115) 
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The probability distributions are then obtained as 

= M _1 [ ] , (116) 

\P3 j 

where 

M = | <>,2 - (^23 + ^2l) <>:<> | . ( 1 17) 

^23 — (S^31 

Using Eq. (|5T|) we can then determine the shear current and the mean shear differences: 

T ( ^i 3 (Ax 3 2 + Ax 2 i) + rii 2 Aa;2i - fii A;E4 3 

n 23 Ax 32 - fi 2 iAx 21 ] . (118) 

fi 34 Ax4 3 - il 32 Ax 32 - ^ 3 iAx 3 i 

By solving Eqs. ()113|) and (|1 14|) numerically for Aq 2 i and Ag 32 given a set of equilibrium rates and distances 
Axji, one can determine the driven transition rates as well as the flux potential for the three state graph of Fig. [9] 
Eqs. (|1 16[) and (|1 18|) can then be used to obtain the probability distributions, the current and the mean shear 
differences. 

The results for the driven transition rates in the forward direction are shown in Fig. 1101 For forward shear three 
of the transition rates are strongly enhanced. However, we observe that the transition 1 — > 2 is attenuated at high 
shear rate, even though it contributes a positive shear increment AX21. This indicates that the system disfavors the 
path via state 2 that requires two transitions to acquire the shear increment A:r 32 . For large shear rates the system 
will thus predominantly choose the high- mobility "channel" 1 — > 3. A similar observation is made for the transition 
2^3, which, for backward driving, remains significant because it connects to the favorable direct channel. 

Accordingly, the steady state occupancies p± and p 3 converge to the value 1/2 for strong driving in both forward 
and backward direction (cf. Fig. [TT]) . On the other hand, the occupancy of state 2, energetically favored over states 
1 and 3 at equilibrium, decays to zero for larger shear. 




IX. SUMMARY 



The search for fundamental principles governing the behavior of systems in out of equilibrium situations has long 
been an area of intensive research. We have discussed a statistical theory (NCDB) that starts from first principles 
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FIG. 11: The steady state distributions pi, p 2 , P3 of Eq. (|116[) plotted parametrically as functions of the current J. Parameter 
values: AE 12 = 3.0, AE 32 = 2.0, wio = 1.0, uis = 1.0, co 32 = 1.5, Ax 21 = 1.8, Ai 32 = 1.0, A143 = 0.5. 



and governs the steady-state motion of any flowing system on which work is done by a weakly- coupled nonequilibrium 
reservoir that is ergodic and microscopically reversible, such as a complex fluid under continuous shear. NCDB 
provides a description of the microscopic dynamics of these systems in the form of exact constraints on the transition 
rates in the driven steady state, arbitrarily far away from equilibrium, akin to the principle of detailed balance for 
equilibrium systems. In this article we have investigated this theory for systems evolving in discrete state spaces. We 
derived a simple relationship between the flux potential and the total exit rates, which leads to a number of important 
further results. 

Two simple sets of invariant quantities have been formulated: (i) the product of forward and reverse transition 
rates and (ii) the differences in total exit rates, for every pair of microstates, equal the corresponding equilibrium 
values and are thus unchanged by the driving. These invariant quantities are non-trivial and experimentally accessible 
predictions of NCDB. We have devised a systematic method to determine the driven transition rates by setting up a 
system of equations from simple network rules. Furthermore, we have investigated properties of the shear current in 
systems satisfying NCDB. We show that the stationary shear current can be determined independently of the steady 
state distribution of microstates. Considering the shear current exhibited by an individual system trajectory, we show 
that the fluctuations of this shear current satisfy a fluctuation relation of the Gallavotti-Cohen type. The validity of 
such a fluctuation relation is related to the Gibbs form of the path distribution p dr (T) (cf. (25j) 

In a system with two types of states all quantities of the NCDB formalism can be expressed analytically. However, 
even for just three states results are only obtainable numerically. Both types of system show consistent behavior. 
In state-spaces containing a loop, the particular path favoured at highest driving is the one most accommodating 
for carrying flux. Even in these simple models the non-local nature of NCDB becomes evident. Whereas mean-field 
theories without the quantity Aq would simply boost any transition in the forward flux direction, NCDB takes into 
account the future propensity for achieving flux. That propensity depends on the global structure of the state space, 
and is communicated in the noise from the reservoir of other systems exploring the possibilities of the steady-state 
dynamics. A particular transition will thus be attenuated if it connects to a state which is blocked or from where 
subsequent transitions carry low flux. This striking property indicates that NCDB might ultimately be able to describe 
the counter-intuitive phase behavior exhibited for example by real complex fluids under shear. 

Future work in this context will predominantly focus on further applications of NCDB to more realistic models as 
well as independent experimental tests. In particular the invariant quantities provide a straightforward criterion 
to check the validity of the theory (cf. [12|]). Even though these new exact relations are based on a rigorous 
statistical mechanical derivation, only comparison with experimental data can shed light on their significance for 
our understanding of nonequilibrium phenomena. 
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